pacman::p_load(sf, spNetwork, tmap, tidyverse)In-class Exercise 3
3 Network Constrained Spatial Point Patterns Analysis
3.1 Loading R packages
3.2 Data Import and Preparation
The code chunk below uses st_read() of sf package to import Punggol_St and Punggol_CC geospatial data sets into RStudio as sf data frames.
#Punggol_St is in ESRI Shapefile format
network <- st_read(dsn = "data/geospatial",
layer = "Punggol_St")Reading layer `Punggol_St' from data source
`C:\zjho008\ISSS626-GAA\In-class_Ex\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TMData has to be in linestring and not multiple linestring. Use st to convert it from multi line to single line.
childcare <- st_read(dsn = "data/geospatial",
layer = "Punggol_CC") %>%
st_zm(drop = TRUE,
what = "ZM") # to remove z valueReading layer `Punggol_CC' from data source
`C:\zjho008\ISSS626-GAA\In-class_Ex\In-class_Ex03\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range: zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TMWe can examine the structure of the output simple data features data tables in R studio.
childcare File is in kml format hence the dimension shown is XYZ (additional dimension). It is important to note the dimension. Upon further inspection under geometry the childcare data has point Z.
Note: for take home exercise under entire data file, 1 folder as rawdata, with another separate folder as data for analysis.
3.3 Visualising the Geospatial Data
Before we jump into the analysis, it is a good practice to visualise the geospatial data. There are at least two ways to visualise the geospatial data. One way is by using plot() of Base R as shown in the code chunk below.
plot(st_geometry(network)) # plotting the road network first, especially when in sf layer
plot(childcare, add = T, col = 'red', pch = 19) # followed by the childcare # since mapped with colours when plotted multiple colours do not appear
# add = T -> T = TRUE the point is plotted twice.Code chunk result when removing the st_geometry:
plot(network)
plot(childcare, add = T, col = 'red', pch = 19)
Note
network has 3 columns: Link ID St_name Geometry
Removing st_geometry will result in individual columns which are pulled out and plotted individually.
To visualise the data with high cartographic quality and in an interactive manner. The mapping function of tmap package can be used as shown in the code chunk below.
tmap_mode('plot')tmap mode set to plottingtm_shape(childcare) + # specifying the layer that is being used
tm_dots(col = "red") +
tm_shape(network) + # to use the extent of the map layer
tm_lines()
tmap_mode('plot')tmap mode set to plottingways to add markers;
https://r-tmap.github.io/tmap/reference/index.html
Specify the shape object: tm_symbols() tm_squares() tm_bubbles() tm_dots() - to keep the size constant when performing zoom functions tm_markers()
Making the plot an interactive layer
tmap_mode('view') #just by switching to 'view' to achieve the interactivitytmap mode set to interactive viewingtm_shape(childcare) + # specifying the layer that is being used
tm_dots(col = "red") +
tm_shape(network) + # to use the extent of the map layer
tm_lines()tmap_mode('plot') # to ensure after the session is ended it will end in the plot mode to reduce resource consumptiontmap mode set to plotting
Note
Childcare & network can be switched on and off accordingly.
3 different data consumptions: 2 layers of ESRI map data (WorldGray Canvas & OpenStreetMap) Topographic layer
mpabox ~ leaflet
While using tmap methods requires a longer code chunk the benefit it brings are the flexibility and customisation that can be done.
3.3.1 Preparing the lixels objects
Before computing NKDE, the Spatial Lines object needs to be cut into lixels with a specified minimal distance. This task can be performed by using lixelize_lines() of spNetwork as shown in the code chunk below.
lixels <- lixelize_lines(network,
700,
mindist = 350)given that it is a road network and in the context of the childcare - so using the reasonable walking distance based on weather and perceived hindrance is about 700 metres based on a study for perceivable walking distance.
mindist is set as half for the minimum walking distance.
2642 segments in the line network. to split into line segment each should be 700 and in the centre the minimum distance should be 350m.
After running the code chunk the segments, the remaining is slightly greater than 350.
if increase to 500 the segment is 2645
if reduce to 150m the segment is still the same 2645.
For take home ex 3, the BMR has to be plotted - a rough gauge of the general distance so we should not use a distance smaller than the point. Calculating the nearest neighbour to find out the nearest neighbour -
based on distances starting on the lowest 25 percentile of accidents along the road segment. Want to acheive a segment that can pick up some accident occurences.
3.3.2 Generating line centre points
Next, lines_center() of spNetwork will be used to generate a SpatialPointsDataFrame (i.e. samples) with line centre points as shown in the code chunk below.
samples <- lines_center(lixels) # sf format3.4 Visualising the lixel segment
tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels) +
tm_lines() +
tm_shape(samples) +
tm_dots(size = 0.01)tmap_mode('plot')tmap mode set to plotting3.4.1 Performing NKDE
We are ready to compute the NKDE by using the code chunk below:
densities <- nkde(network,
events = childcare,
w = rep(1, nrow(childcare)),
samples = samples,
kernel_name = "quartic",
bw = 300,
div = "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE)# avoid gaussian if intensity changes to negative # 3 methods: simple, continous, discontinous
the computed density values (i.e. densities) into samples and lixels objects as density field.
samples$density <- densities
lixels$density <- densitiesTo append the intensity values into the simple tibular frame or lixel data frame simialr to a left join.
Avoid sorting to avoid changing the sequence.
values attached to the line and point.
Since the svy21 projection system is in metres, the computed density values are very small i.e. 0.0000005. The code chunk below is used to rescale the density values from number of events per metre to number of events per kilometre.
# rescaling
samples$density <- samples$density*1000
lixels$density <- lixels$density*1000The code below uses appropriate functions of tmap package to prepare interactive and high cartographic quality map visualisation.
tmap_mode(mode = c("view"))tmap mode set to interactive viewingtm_shape(lixels) +
tm_lines(col = "density") +
tm_shape(childcare) +
tm_dots()tmap_mode("plot")tmap mode set to plottingkfun_childcare <- kfunctions(network,
childcare,
start = 0,
end = 1000,
step = 50,
width = 50,
nsim = 99, # simulations are starting from zero
resolution = 50,
verbose = FALSE,
conf_int = 0.05)The output of kfunctions() is a list with the following values:
- plotkA, a ggplot2 object representing the values of the k-function
- plotgA, a ggplot2 object representing the values of the g-function
- valuesA, a DataFrame with the values used to build the plots
Visualising the ggplot2 object of k-function by using the code chunk below.
::: panel-tabset ## K-Function
kfun_childcare$plotk # whether to plot G / K function
2 possible patterns observed
regular pattern below the envelope - showing signs of regularity - childcare centres near to each other e.g. at 200m apart which is showing the signs of regularity
and complete spatial randomness at the upper portion.
3.5 G-Function
kfun_childcare$plotg # whether to plot G / K function
both functions are returned